df_od <- read_csv(here::here("data", "VSRR_Provisional_Drug_Overdose_Death_Counts.csv"))
df_od$Date <- parse_date_time(paste(tolower(df_od$Month), df_od$Year), orders=c("by"))
df_od$NumDeaths <- df_od$`Data Value`
df_od <- df_od %>% filter(`Indicator`=='Number of Deaths' & `State` != 'US')
df_od <- df_od %>% left_join(., pop, by=c("State"="State"))
df_od$FatalitiesPerCap <- df_od$`Data Value`/df_od$pop*100
p1 <- ggplot(df_od, aes(x=Date, y=`FatalitiesPerCap`, group=`StateName`, colour=`StateName`)) +
geom_line() +
gghighlight(max_highlight = 8L, max(`FatalitiesPerCap`), use_direct_label=FALSE, label_key=`FatalitiesPerCap`) +
labs(title="WEST VIRGINIA'S OPIOID PROBLEM STANDS \nIN STARK CONTRAST TO REST OF COUNTRY",
subtitle = "DRUG OVERDOSES PER CAPITA OVER THE YEARS",
caption="SOURCE: NATIONAL CENTER FOR HEALTH STATISTICS",
x="DATE", y="% FATAL DRUG OVERDOSES PER CAPITA") +
jtheme + jfill + jcolor
p1
Arguably no one has been hit harder by the opioid epidemic than those living in the rust belt and the rural south. One can see above that rural America (Alabama, Arkansas, Kentucky, Tennesee, and West Virginia) comprise some of the states hit hardest per capita by opioid related deaths. Many factors contribute to this scourge on the region, including economic disenfranchisement and high rates of disability enrollment.
df_ods_by_state <- aggregate(df_od, by=list(df_od$StateName), FUN=mean, na.rm=TRUE)
df_ods_by_state$StateName <- df_ods_by_state$Group.1
df_ods <- df_ods_by_state %>% select(StateName, FatalitiesPerCap)
mean_ods <- mean(df_ods$FatalitiesPerCap)
df_ods$Above <- df_ods$FatalitiesPerCap > mean_ods
df_ods$Above[df_ods$Above == TRUE] <- "BLUE"
df_ods$Above[df_ods$Above == FALSE] <- "RED"
p2 <- ggplot(df_ods, aes(x=reorder(StateName, FatalitiesPerCap), y=FatalitiesPerCap, label=format(FatalitiesPerCap, digits=2, nsmall=2))) +
geom_point(stat='identity', aes(col=Above), size=6) +
scale_color_manual(name="FATALITIES BY OVERDOSE (% OF POP)",
labels = c("BELOW AVERAGE", "ABOVE AVERAGE")) +
labs(title="NEW YORK AND WEST VIRGINIA LIE ON OPPOSITE \nPOLES OF THE OPIOID EPIDEMIC",
subtitle="FATAL DRUG OVERDOSES BY STATE",
caption="SOURCE: NATIONAL CENTER FOR HEALTH STATISTICS",
x="STATE", y="FATALITIES PER CAPITA (% OF POP)") +
ylim(0, 2.0) +
coord_flip() +
jtheme + jfill + jcolor
p2
While West Virginia is most intensely affected among US states, no state is immune to the epidemic; even in New York, the state whose opioid-related death rate is lowest per capita, one half of one percent of the population will die from a drug-related overdose. The standard deviation of opioid abuse’s impact on states is smaller than one might guess if they are only reading articles (https://www.nytimes.com/interactive/2018/us/west-virginia-opioids.html) about the plight of West Virginia. While West Virginia is an outlier as shown in the above plot, this is a nationwide epidemic.
DRUG_NAMES <- c("VICODIN", "OXYCODONE", "OXYCONTIN", "PERCOCET",
"OPANA", "KADIAN", "AVINZA", "FENTANYL")
fname_drug_utilization <- "State_Drug_Utilization_Data_"
# First time processing of DataFrame
#df_util <- read_csv(here::here("data", paste(fname_drug_utilization, "2008", ".csv", sep=""))) %>% filter(`Product Name` %in% DRUG_NAMES)
#for (i in 2009:2018) {
# df_util <- rbind(df_util, read_csv(here::here("data", paste(fname_drug_utilization, i, ".csv", sep=""))) %>% filter(`Product Name` %in% DRUG_NAMES))
#}
#readr::write_csv(df_util, here::here("data", "df_util.csv"))
df_util <- read_csv(here::here("data", "df_util.csv"))
df_rx <- aggregate(df_util[,c("Number of Prescriptions", "Total Amount Reimbursed")], by=list(df_util$Year, df_util$`Product Name`), FUN="sum", na.rm=TRUE)
df_rx$`Drug Name` <- df_rx$Group.2
p3 <- df_rx %>% ggplot(aes(Group.1, `Total Amount Reimbursed`, size=`Number of Prescriptions`, colour=`Drug Name`)) +
geom_point(alpha = 0.7) +
geom_line(size=1) +
scale_y_continuous(labels = dollar) +
labs(title = "PAINKILLER REIMBURSEMENTS SHRINK WHILE \nPRESCRIPTIONS REMAIN FAIRLY CONSTANT",
subtitle = "EXPENSE REIMBURSEMENT & PRESCRIPTIONS OVER TIME",
caption = "SOURCE: DATA.MEDICAID.GOV",
x = 'YEAR', y = 'TOTAL $ REIMBURSED') +
jtheme + jfill + jcolor
p3$labels$size <- "NUMBER OF PRESCRIPTIONS"
p3$labels$colour <- "DRUG NAME"
p3
2015 was a clear reckoning for the US Government’s reimbursement of presciption opioids, those contributing most towards the opioid epidemic. The data above shows the prescriptions in cicle size over time, along with the total amount reimbursed, thanks to medicaid data. Economists might anticipate that as less and less reimbusement money is available for prescription opioids, doctors would prescribe fewer pills. Though we do see slightly larger circles as reimbursements go up, it seems that prescriptions aren’t very responsive to drug reimbursement. This relationship may show different results in comparing reimbursement to the number of filled prescriptions at the pharmacy.
df_util$StateName <- abbr2state(df_util$State)
# Drug Overdoses
DATA_OVERDOSE <- "VSRR_Provisional_Drug_Overdose_Death_Counts.csv"
df_od <- read_csv(here::here("data", DATA_OVERDOSE))
df_od$Date <- parse_date_time(paste(tolower(df_od$Month), df_od$Year), orders=c("by"))
df_od$NumDeaths <- df_od$`Data Value`
df_od <- df_od %>% filter(`Indicator`=='Number of Deaths' & `State` != 'US')
df_od <- df_od %>% left_join(., pop, by=c("State"="State")) %>% left_join(., df_util, by=c("State"="State", "Year"="Year"))
df_od$FatalitiesPer <- df_od$`Data Value`*1000/df_od$pop
df_od$`Number of Prescriptions` <- df_od$`Number of Prescriptions`*1000/df_od$pop
yr=2018
df_ods_by_state <- df_od %>% filter(Year == yr)
df_ods_by_state <- aggregate(list(FatalitiesPer=df_ods_by_state$FatalitiesPer, NumRxPer=df_ods_by_state$`Number of Prescriptions`), by=list(StateName=df_ods_by_state$StateName.x), FUN="sum", na.rm=TRUE)
# Import Cartogram Polygons, merge drug overdose data
spdf <- geojson_read(here::here("data", "us_states_hexgrid.geojson"), what = "sp")
spdf@data = spdf@data %>% mutate(google_name = gsub(" \\(United States\\)", "", google_name))
spdf@data = spdf@data %>% left_join(., df_ods_by_state, by=c("google_name"="StateName"))
cartogram <- cartogram_cont(spdf, 'FatalitiesPer')
carto_fortified <- tidy(cartogram, region = "google_name") %>% left_join(. , cartogram@data, by=c("id"="google_name"))
centers <- cbind.data.frame(data.frame(gCentroid(cartogram, byid=TRUE), id=cartogram@data$iso3166_2))
p_green <- ggplot() +
geom_polygon(data = carto_fortified, aes(fill = FatalitiesPer, x = long, y = lat, group = group) , size=0.05, alpha=0.9, color="#f5f5f9") +
scale_fill_gradientn(colours=brewer.pal(9, "Greens"), name="OPIOID RELATED DEATHS", guide=guide_legend( keyheight = unit(3, units = "mm"), keywidth=unit(15, units = "mm"), title.position = 'top', label.position = "bottom") ) +
geom_text(data=centers, aes(x=x, y=y, label=id), color="#444444", size=3, alpha=0.6) +
theme_void() +
labs( title=paste("WHERE OVERDOSE DEATHS HIT HARDEST\n",yr),
subtitle="CARTOGRAM OF DEATHS ACROSS THE US",
caption="SOURCE: VSSR, DATA.MEDICAID.GOV") +
jtheme_nogrid +
coord_map()
p_blue <- ggplot() +
geom_polygon(data = carto_fortified, aes(fill = NumRxPer, x = long, y = lat, group = group) , size=0.05, alpha=0.9, color="#f5f5f9") +
scale_fill_gradientn(colours=brewer.pal(9, "Blues"), name="PRESCRIPTION RATES", guide=guide_legend( keyheight = unit(3, units = "mm"), keywidth=unit(15, units = "mm"), title.position = 'top', label.position = "bottom") ) +
geom_text(data=centers, aes(x=x, y=y, label=id), color="#444444", size=3, alpha=0.6) +
theme_void() +
labs( title=paste("WHERE DRUGS LEAVE THE HOSPITAL\n",yr),
subtitle="CARTOGRAM OF KEY OPIOID PRESCRIPTION RATES ACROSS THE US",
caption="SOURCE: VSSR, DATA.MEDICAID.GOV") +
jtheme_nogrid +
coord_map()
#grid.arrange(p_green, p_blue, nrow = 2)
p_green
The rust belt in 2018 saw the worst of opioid related drug overdoses. This likely has to do with smuggling routes, with major hubs in cities like Chicago. The rural west appears to be relatively minimally affected, with states along the Sierra Nevadas showing low rates of overdose deaths. In the Northeast, results are mixed, with coastal mid-atlantic states faring much better than the rural border states like New Hampshire.
p_blue
Correlation to states with high prescription rates doesn’t tell much of a story when comparing cartograms. Maryland had by far the most opioid prescriptions, followed by a cluster of states in the Southwest, omitting Utah, whose large Mormon population dampens its statistics related to drug utilization. There doesn’t appear to be a relationship within given years between opioid prescription and abuse within states. This could hint at a black market, or that most people are not getting addicted to opioids from injury.
US <- map_data("world") %>% filter(region=="USA")
data=us.cities
df_util$`NUMBER OF PRESCRIPTIONS` = df_util$`Number of Prescriptions`
df_util$`PRODUCT NAME` = df_util$`Product Name`
unique(df_util$`Product Name`)
## [1] "KADIAN" "AVINZA" "OXYCODONE" "OXYCONTIN" "PERCOCET" "FENTANYL"
## [7] "OPANA" "VICODIN"
df_util %>% filter(df_util$`Product Name` %in% c("OXYCODONE", "OXYCONTIN")) %>%
#arrange(`Number of Prescriptions`) %>%
#mutate( name=factor(name, unique(name))) %>%
#mutate(pop=pop/1000000) %>%
ggplot() +
geom_polygon(data = US, aes(x=long, y = lat, group = group), fill="grey", alpha=0.3) +
geom_point(width=2, height=2, aes(x=Longitude, y=Latitude, size=`NUMBER OF PRESCRIPTIONS`, color=`PRODUCT NAME`), shape=20, stroke=FALSE) +
scale_size_continuous(name="NUMBER OF PRESCRIPTIONS", range=c(1,12)) +
#scale_alpha_continuous(name="Year", range=c(0.1, .7)) +
xlim(-125, -70) + ylim(25,50) + coord_map() +
guides(color=guide_legend(title="PRODUCT NAME")) +
jtheme_nogrid + jcolor +
labs(title = "OXY'S MOVEMENT SHOWS A LULL\nBETWEEN 2012-2014", subtitle='PRESCRIPTIONS THROUGH THE YEARS: {frame_time}', caption="SOURCE: MEDICAID.GOV") +
#transition_states(year, transition_length = 2, state_length=3) +
transition_time(Year) +
ease_aes('linear', interval = 1.0)